Classifying the expansion kinetics and critical surface dynamics of growing cell 

populations 
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Based on a cellular automaton model the growth kinetics and the critical surface dynamics of 
cell monolayers is systematically studied by variation of the cell migration activity, the size of the 
proliferation zone and the cell cycle time distribution over wide ranges. The model design avoids 
lattice artifacts and ensures high performance. The monolayer expansion velocity derived from our 
simulations can be interpreted as a generalization of the velocity relationship for a traveling front 
in the Fisher-Kolmogorov-Petrovskii-Piskounov (FKPP) equation that is frequently used to model 
tumor growth phenomena by continuum models. The critical surface dynamics corresponds to the 
Kardar-Parisi-Zhang (KPZ) universality class for all parameters and model variations studied. While 
the velocity agrees quantitatively with experimental observations by Bru et al, the critical surface 
dynamics is in contrast to their interpretation as generic molecular-beam-epitaxy-like growth. 
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Model simulations of tumor growth and therapy have 
attracted wide interest An important issue to 

which models can contribute is the classification of the 
tumor growth pattern by generic mechanisms at the level 
of the individual cell actions (migration, division etc.). 
These actions subsume the effect of the molecular inter- 
and intra-cellular regulation. The models can serve to 
identify those cell activities that would result in a maxi- 
mal inhibition of multi-cellular growth and invasion, and 
thereby help to identify possible molecular drug targets. 
Bru et al 4] analyzed the growth kinetics and critical 
surface dynamics of many different tumors in-vitro and 
in-vivo. They quantified the dynamics of the tumor sur- 
face by three critical exponents used to classify crystal 
growth phenomena into universality classes |5j. They 
found a generic linear growth phase of in-vitro growing 
cell lines for large cell populations and a molecular-beam- 
epitaxy (MBE)-like dynamics of the tumor surface both 
in-vitro and in-vivo. They proposed a tumor therapy 
based on these findings |6(. 

In this letter we analyze a class of cellular automaton 
(CA) tumor growth models on an irregular lattice by ex- 
tensive computer simulations. CA tumor growth models 
enjoy wide interest since they permit to represent each 
cell individually at moderate computational expense. In 
our model cells can divide, push neighbor cells and mi- 
grate. The choice of the model rules is guided by com- 
parison with an off-lattice model. By using the irregular 
lattice we ensure isotropy and homogeneity of space, and 
cell sizes that are sharply peaked around a prescribed 
average value. Both the expansion speed and the spatial 
pattern formed differ from results on a periodic lattice. 
We systematically analyze our growth model with respect 



to the hopping rate, proliferation depth and dispersion of 
the cell cycle time distribution and show that the expan- 
sion dynamics can be mapped onto the functional form 
of the traveling wave velocity of the Fisher-Kolmogorov- 
Petrovskii-Piskounov (FKPP) equation [7|. The model 
reproduces the monolayer expansion kinetics experimen- 
tally found by Bru Q. The critical surface growth dy- 
namics suggests a Kadar-Parisi- Zhang (KPZ)-like Q be- 
havior over a wide range of parameters and for varying 
cell migration mechanisms, supporting the critical com- 
ment by Buceta and Galeano [|| on the conjecture by 
Bru et. al. ;4J. Our findings comply with the results in 
the classical Eden model [T(J- 
Our model is based upon the following assumptions: 
[Rl] Lattice generation: Starting from a regular square 
lattice with spacing I, an irregular lattice 7\j is gener- 
ated by Delauney triangulation. A biological cell is rep- 
resented as shown in Fig^a) (white). 
[R2] Exclusion principle: Each lattice site can be occu- 
pied by at most one single cell. 

[R3] Cycle time: The cell cycle time t' is Erlang dis- 
tributed (with a parameter m): 



Z(t') = A, 
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with A m = m such that (r') = r = 1. 
[i24] Proliferation depth: A cell can divide if and only if 
there is at least one free neighbor site within a circle of 
radius AL around the dividing cell (Fig. 2] (a) , green) . 
[RS\ Cell migration: We consider three alternative migra- 
tion rules: R5(i) A cell moves with rate <fi to a free neigh- 
bor site, irrespectively of the number of neighbor cells 
before and after its move. This rule corresponds to the 
case of no cell-cell adhesion. R5(ii) Cells move with rate 4> 
if by this move the cell is not isolated. R5(iii) Cells move 
with a rate ob exp{~ AE / F T } with AE = E(t+At)-E(t), 
where At is the time step, E(t) is the total interaction 
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Figure 1: (Color online) (a) Construction of the CA lattice: 
One point (black, green) is placed in every square of a square 
lattice at a random position r i . A Voronoi tesselation is con- 
structed form these points such that each cell consists of all 
points in space that are closer to the lattice point than to 
any other r_ k . The shape of a biological cell (white) is iden- 
tified with the corresponding Voronoi polygon (blue lines). 
Polygons that share a common edge are defined as neighbor- 
ing and connected by red lines (Delauney triangulation) . (b) 
Probability density distribution of the cell area for the CA- 
lattice in (a) (brown) and for a random initial distribution 
of points (red), (c) Cell cluster morphology for m — 10 4 , 
AL — 1 on (i) the CA lattice in (a), (ii) square-, (iii) hexag- 
onal lattice, (iv) lattice with Moore neighborhood (nearest 
neighbors along the axes and the diagonals), (v) off-lattice 
cluster HE3. 

energy of the multi-cellular configuration, Ft ~ 10 -16 J 
is a "metabolic" energy £2, AE/F T ~ 0(1)- O(10) 0. 
This induces migration towards locations with a larger 
number of neighbor cells. 

By [Rl] we generate an unstructured lattice with a sym- 
metric cell area distribution sharply peaked around its 
average A = I 2 (see Fig|U(a),(b)). [R3] considers that 
experiments indicate a T-like distribution of the cell cy- 
cle controlled by cell cycle check points |l2j. [R4] takes 
into regard that the growth speed of tumors is usually 
incompatible with the assumption that only cells at the 
border are able to divide (as in the Eden model see 
|3). Therefore we assume that a dividing cell is able to 
trigger the migration of at most k neighbor cells into the 
direction of minimum mechanical stress (see Fig^ (a)). 
If a cell divides, one of its daughter cells is placed at 
the original position, the other cell is placed next to it 
and the local cell configuration is shifted and re-arranged 
along the line that connects the dividing cell with the 
closed free lattice site within a circle of radius AL such 
that the latter is now occupied (see FigU(a)). This al- 
gorithm mimics a realistic re-arrangement process that 
may occur from active cell migration as a response to a 
mechanical stimulus, cf. Ref. |14j. Isolated cells perform 
a random- walk-like motion (e.g. |l5j). We consider dif- 
ferent migration rules R5(i)-(iii) to comprise a class of 
potential models with biologically realistic behavior. 
The model parameters are the average cell cycle time r 



and its distribution /(r') controlled by the parameter to, 
the migration rate </>, the proliferation depth AL, and, in 
case of an energy-activated migration rule, the e nerg y E. 
Programmed cell death can easily be integrated [21j but 
is omitted here. Rules [R1-R5] can be formalized by the 
master equation 

d t p(Z,t)= W z >->zp{Z',t) -W z ^z-p{Z,t). (2) 

Z'^Z 

Here p(Z, t) denotes the multivariate probability to find 
the cells in configuration Z and W{Z' — > Z) denotes the 
transition rate from configuration Z' to configuration Z. 
A configuration Z = {...,Xi—\,Xi,x%+\, ...} consists of lo- 
cal variables xi = {0, 1} with Xi = if lattice site i is 
empty, and X\ = 1 if it is occupied by a cell. For the 
simulation we use the Gillespie algorithm |l6j . i.e, the 
time-step of the event-based simulation is a random num- 
ber given by At = —yk-ln(l — £). Here, £ is a random 
number cquidistributcd in [0,1), Wz = J^Z'^z'^z is 
the sum of all possible events which may occur at time t. 
Here we assume that the rate at which a cell changes its 
state by a hop, a progress in the cell cycle, or a division 
is independent of the number of accessible states as long 
as at least one state, that is, one free adjacent lattice site 
in case of a hop and one free site within a circle of radius 
AL in case of a division, is accessible. This may be justi- 
fied by noting that cells - in contrast to physical particles 
- are able to sense their environment and therefore the 
direction into which they can move. 
We analyze the growth kinetics by the cell population size 
N(i) (number of cells at time t) and the radius of gyration 

Rgyr(t) = \f jfJ2?=i(Lt- R a ) 2 - Here R = ± J^Zi U is 
the position of the center of mass. For a compact circular 
cell aggregate (in d = 2 dimensions), R gyr is related to 
the mean radius R(t) = ^ J 2 ™ R(<p, t)d(f (polar angle <p) 
of the aggregate by R — R gyr \/2. 

To interpret the rules and parameters of the CA model 
in terms of growth mechanisms we compare it with the 
stochastic single-cell-based off-lattice growth model in 
Ref. (Fig. |2l . In this model cell motion contains 
an active random component and a component triggered 
by mechanical forces between cells, and between cells and 
the substrate 01 ■ During cell division the cell gradually 
deforms and divides into two daughter cells as long as the 
degree of deformation and compression is not too large. 
As illustrated in Fig. [2] the lattice model is able to cap- 
ture the behavior of the off-lattice model and agrees with 
the experimental findings in Refs. Q provided the pa- 
rameters AL, 4>, r, to are chosen properly. AL controls 
the effective thickness of the proliferative rim; in the off- 
lattice model it depends on the mechanisms that control 
the proliferation by contact inhibition, on the material 
properties of the cell (the Young modulus, the Poisson 
number etc.), and on the ability of a cells to move in re- 
sponse to a mechanical stimulus 3] . 
At large m the tumor border becomes smoother and the 
tumor shape reflects the symmetry of the underlying lat- 
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Figure 2: (Color online) (a) Mean radius R of the cell aggre- 
gate vs. time t. Full circles: experimental findings for C6 rat 
astrocyte glioma cells (0). (b) Cell cycle time distribution 
/(r') for the off-lattice model and the CA growth model in 
comparison with the Erlang distribution, (m = 60, AL — 9, 
= 0) 
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tice (Fi g. [D (c)(ii-iv)); this effect is known as noise re- 
duction |l9j • Such lattice-induced asymmetries could sig- 
nificantly disturb the analysis of the surface growth dy- 
namics in circular geometries. We have chosen a Voronoi 
tesselation, in which such artifacts do not occur (Fig. 
H (a) ,(c) (i)) . Fig. shows a systematic study of the 
growth kinetics for free hopping (Rule R5(i)). All quan- 
tities are plotted in multiples of r and I, which are the 
reference time and length scale, respectively. Initially, 
the cell population size grows exponentially fast with 
NU) = N(0)exp(t/T cS ) where r^ 1 = (2 1 /" 1 - ljmr 1 
|l8| . The duration of the initial phase increases with 
AL and <fi. The growth law for the diameter depends 
on (/). If <f> = 0, the initial expansion of the diameter 
is exponentially fast, too. If <fr > 0, cells initially de- 
tach from the main cluster and the diameter grows dif- 
fusively, with L = 2\f2R gyr oc \[M$> + Tpr^gji where 
A sa 1.2 is a lattice-dependent fit constant (Fig. Ela)). 
For t/r < 2, R gyr oc t (Fig. OJa)). This regime disap- 
pears for N(0) ^> 1 (see [l8|). As soon as cells in the 
interior of the aggregate are incapable of further division 
the exponential growth crosses over to a linear expansion 
phase. Fig. 01 shows v 2 vs. (b) (AL) 2 , (c) (f>, and (d) m 
for large N (N ~ 10 5 cells). The model can explain the 
experimentally observed velocity-range in Ref. Q. As 
t — ► oo, L = v(m, <f>, AL)t with 



B 2 ([AL'(AL)]Vr l W/T eff ), 



(3) 



B ps 1.4 (lines in Fig. EJd-c). AL'(AL) (w 1+0.6(AL-1)) 
results from the average over all permutations to pick 
boundary cells within a layer of thickness AL. For 
AL/t c s -c VW^ff e 1 n - El has the same form as for 
the FKPP equation, (e.g. 0]). 

Next, to determine the universality class we deter- 
mine the roughness exponent a and the dynamic expo- 
nent z from the dynamic structure function S(k,t) = 
(R(k, t)R(—k, t)) where R(k, t) is the Fourier transform 
of the local radius R(s,t) and (...) denotes the aver- 
age over different realizations of the growth process (e.g. 
[23). Here s is the arclength as in Ref. 0. The third 
exponent, the growth exponent (3, can be obtained from 



Figure 3: (Color online) (a) Y = Rg yr /((f> + 1/t c b) vs. t/r 
for m — 0, AL — 1 and different values for <f>. (b-d): Growth 
in the linear expansion regime (N ~ 10 s ). (b) Square of 
expansion velocity, v 2 , vs. square of the proliferation zone, 
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m — 0). (c) v 2 vs. <j> (triangles: AL — 1, circles: AL — 3, 
squares: AL — 6, stars: AL = 10; m = 0). (d) v vs. m 
(AL = 1, (j> = 0). The lines are fits using eqn. 




Figure 4: (Color online) (a) Dynamic structure function 
S(k,t) vs. k for [R5(i)], AL = 0, <f> = 0, m = 0. Inset: 
rescaled function S(k,t)k 2a+1 vs. kt 1/z (a = 0.5, z = 3/2). 
(b) S(k, t) vs. k for four alternative parameter sets: (A) tri- 
angles: m = 5 (AL = 0, <t> = 0), (B) circles: AL = 6 (m = 0, 
= 0), (C) stars: R5(ii) = 100 (m = 0, AL = 0), (D) 
squares: R5(iii) AE/Ft = Eo + n ■ Eb with Eb = 10, n 
neighbors, Eo = 5 surface binding (m = 0, AL — 0) 20]. The 
dashed lines are guides to the eye showing a — 0.5. 
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the scaling relation (3 = a/z. In test simulations compar- 
ing constant angle segments Aip with constant arclength 
intervals As we did not find noteworthy differences. For 
sclf-affme surfaces in absence of any critical length-scale 
the dynamic structure function has the Family- Vicsek 
scaling form 

S(k,t) = fc-( 2Q+1) s(H 1 / 2 ) (4) 

I / const - if u > 1 

s{u = kt') = | u _ (aa+1) f/ u<<1 (5) 

At u = 1 a crossover occurs. For «»1 curves measured 
at different times collapse onto a single line; at u <C 1 
they split. We have calculated S(k,t) for rules R5(i) 
and <j) > 0, R5(ii) and RS(iii) (Fig. Q). The final cell 
population size was of O(10 5 ) cells which is the typical 
size of the cell populations in Ref. All these results 
suggest KPZ-like dynamics with a — 1/2, z — 3/2 and 
(3 = 1/3 rather than the MBE universality class, i.e., 
critical exponents a = 3/2, z = 4 and /3 = 3/8 inferred 
in 0. The parameter range of <fr G [0,100) captures 
most cell lines studied in Ref. £| (for I = 10/zm, 
r = 24ft,, cf) = 100 corresponds to a diffusion constant of 
D = 10- 10 cm 2 /s). 

In conclusion we have analyzed the expansion ki- 
netics and critical surface dynamics of two-dimensional 
cell aggregates by extensive computer simulations within 
a CA model which avoids artifacts from the symmetry 



of regular lattices. The growth scenarios are compatible 
with experimental observations. The asymptotic expan- 
sion velocity has a form that is reminiscent of the front 
velocity of the FKPP equation. The same expansion 
velocity can be obtained for different combinations of the 
migration and division activities of the cell and of the 
cycle time distribution. Recently, mathematical models 
based on the FKPP equation were used to predict the 
distribution of tumor cells for high-grade glioma in re- 
gions which are below the detection threshold of medical 
image techniques [24j. We believe such predictions must 
fail since the FKPP equation lacks some important 
parameters such as the proliferation depth which is 
why it is not sensitive to relative contributions of the 
proliferation depth and free migration. We observed in 
our simulations that these relative contributions in fact 
determine the cell density profile at the tumor-medium 
interface: the larger the fraction of free migration is, the 
wider is the front profile even if the average expansion 
velocity is constant. 

The critical surface dynamics found in our simulations 
does not comply with the interpretation of experimental 
observations by Bru et. al. 0] even for the migration 
mechanism they suggested (R5(iii)). We propose to 
re-analyze the corresponding experiments and track the 
paths of marked cells. 
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